Most data has a spatial dimension to it. Knowing ‘where’ the data is coming from is often as crucial as knowing the ‘what’, ‘when’ and ‘who’ dimensions of a given dataset. Therefore it should be no surprise that R has a rich suite of packages for constructing maps and analyzing spatial data. R’s capability has grown so much over the years that it’s functionality rivals many dedicated geographic information systems (GIS). During this Meetup the basics for managing and mapping spatial data with be introduced, using the following packages: sf, ggplot2, tmap, mapview and leaflet.
How you plot your data in R depends on what format it’s in. R has several different formats for managing spatial data (e.g. sp vs sf), and different formats are only compatible with certain plotting systems. Thankfully converting between different spatial formats is not difficult.
library(maps)
map(database = "state")
map(database = "state", region = "indiana")
library(USAboundaries)
library(sf)
cnty <- us_counties()
cnty <- subset(cnty, !state_name %in% c("Alaska", "Hawaii"))
vars <- c("statefp", "countyfp")
plot(cnty[vars])
https://github.com/potterzot/rnassqs
# devtools::install_github('potterzot/rnassqs')
# library(rnassqs)
# source("C:/Users/steph/Nextcloud/code/api_keys.R")
data(state)
st <- state.abb
# corn_us <- lapply(st, function(x) {
# cat("getting", x, as.character(Sys.time()), "\n")
# tryCatch({
# corn = nassqs_yield(
# list("commodity_desc"="CORN",
# "agg_level_desc"="COUNTY",
# "state_alpha"=x
# ),
# key = nass_key
# )},
# error = function(err) {
# print(paste("Error occured: ",err))
# return(NULL)
# }
# )
# })
# corn_us <- do.call("rbind", corn_us)
#
# save(corn_us, file = "C:/workspace2/corn_us.RData")
# write.csv(corn_us, file = "nass_corn_us.csv", row.names = FALSE)
load(file = "C:/Users/Stephen.Roecker/Nextcloud/data/corn_us.RData")
corn_yield <- subset(corn_us, short_desc == "CORN, GRAIN - YIELD, MEASURED IN BU / ACRE")
corn_yield <- within(corn_yield, {
Value = as.numeric(Value)
year = as.numeric(year)
state_name = NULL
state = state_alpha
})
cnty_corn <- merge(cnty, corn_yield,
by.x = c("state_abbr", "countyfp"),
by.y = c("state_alpha", "county_code"),
all.x = TRUE
)
corn_states <- c("IL", "IA", "IN", "MI", "MN", "MO", "NE", "OH", "SD", "ND", "WI")
library(dplyr)
library(ggplot2)
group_by(corn_yield, state_alpha, year) %>%
summarize(
yield_low = min(Value, na.rm = TRUE),
yield_median = median(Value, na.rm = TRUE),
yield_max = max(Value, na.rm = TRUE)
) %>%
filter(state_alpha %in% corn_states) %>%
mutate(state = state_alpha,
source = "NASS"
) %>%
ggplot() +
geom_line(aes(x = year, y = yield_median, col = source)) +
geom_ribbon(aes(x = year, ymin = yield_low, ymax = yield_max), alpha = 0.25) +
facet_wrap(~ state) +
ylab("yield per county (bu/acre)") +
ggtitle("USDA-NASS Corn Yields")
# geom_point(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, col = "NASIS"), size = 1) +
# geom_ribbon(data = yld_sum2[yld_sum2$state %in% corn_states, ], aes(x = year, ymin = yield_low2, ymax = yield_max2, col = "NASIS"), alpha = 0.25) +
# geom_pointrange(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, ymin = yield_low2, ymax = yield_max2, col = source))
IN <- sf::read_sf(dsn = "D:/geodata/soils/soils_GSMCLIP_mbr_2599033_03/wss_gsmsoil_IN_[2006-07-06]/spatial/gsmsoilmu_a_in.shp", layer = "gsmsoilmu_a_in")
# simplify polygons
IN <- rmapshaper::ms_simplify(IN)
library(aqp)
library(soilDB)
# download lab locations for the Miami soil series
miami <- fetchKSSL("Miami")
miami <- site(miami)
miami <- subset(miami, complete.cases(x, y))
miami <- within(miami, {
lon = x
lat = y
})
str(miami)
## 'data.frame': 56 obs. of 34 variables:
## $ pedon_key : chr "7153" "7154" "8836" "8840" ...
## $ pedlabsampnum : chr "79P0060" "79P0061" "81P0725" "81P0729" ...
## $ pedon_id : chr "78IN177007" "78IN177008" "S1981IN011002" "S1981IN011006" ...
## $ pscs_top : int 15 30 20 20 25 15 NA NA NA NA ...
## $ pscs_bottom : int 65 80 70 70 75 65 NA NA NA NA ...
## $ pedon_completeness_index: int 1 1 7 7 6 7 6 6 6 6 ...
## $ state : chr "Indiana" "Indiana" "Indiana" "Indiana" ...
## $ mlra : chr "111A" "111A" "111A" "111A" ...
## $ geomposhill : chr "" "" "" "" ...
## $ geomposmntn : logi NA NA NA NA NA NA ...
## $ geompostrce : logi NA NA NA NA NA NA ...
## $ geomposflats : logi NA NA NA NA NA NA ...
## $ hillslopeprof : chr "" "" "" "" ...
## $ geomslopeseg : chr "" "" "" "" ...
## $ pmgroupname : chr "" "" "" "" ...
## $ drainagecl : chr "" "" "Well drained" "Well drained" ...
## $ taxonname : chr "Miami" "Miami" "Miami" "Miami" ...
## $ taxorder : chr "Alfisols" "Alfisols" "Alfisols" "Alfisols" ...
## $ taxsuborder : chr "Udalfs" "Udalfs" "Udalfs" "Udalfs" ...
## $ taxgrtgroup : chr "Hapludalfs" "Hapludalfs" "Hapludalfs" "Hapludalfs" ...
## $ taxsubgrp : chr "Typic Hapludalfs" "Typic Hapludalfs" "Oxyaquic Hapludalfs" "Oxyaquic Hapludalfs" ...
## $ taxpartsize : chr "fine-loamy" "fine-loamy" "fine-loamy" "fine-loamy" ...
## $ taxpartsizemod : chr "" "" "" "" ...
## $ taxceactcl : chr "" "" "active" "active" ...
## $ taxreaction : chr "" "not used" "" "" ...
## $ taxtempcl : chr "mesic" "" "mesic" "mesic" ...
## $ taxmoistscl : chr "" "" "Oxyaquic" "Oxyaquic" ...
## $ taxtempregime : chr "mesic" "mesic" "mesic" "mesic" ...
## $ soiltaxedition : chr "" "" "tenth edition" "tenth edition" ...
## $ osdtypelocflag : int 0 0 0 0 0 0 0 0 0 0 ...
## $ x : num -85 -85 -86.3 -86.3 -86.3 ...
## $ y : num 39.8 39.8 40 40 40 ...
## $ lat : num 39.8 39.8 40 40 40 ...
## $ lon : num -85 -85 -86.3 -86.3 -86.3 ...
# construct
# sp object
library(sp)
miami_sp <- SpatialPointsDataFrame(
data = miami,
coords = cbind(miami$lon, miami$lat),
proj4string = CRS("+init=epsg:4326")
)
# sf object
library(sf)
miami_sf <- st_as_sf(
miami,
coords = c("lon", "lat"),
crs = 4326
)
# data structures
str(miami_sp, 2)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 56 obs. of 34 variables:
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:56, 1:2] -85 -85 -86.3 -86.3 -86.3 ...
## .. ..- attr(*, "dimnames")=List of 2
## ..@ bbox : num [1:2, 1:2] -89.3 39 -83.9 42.4
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
str(miami_sf)
## Classes 'sf' and 'data.frame': 56 obs. of 33 variables:
## $ pedon_key : chr "7153" "7154" "8836" "8840" ...
## $ pedlabsampnum : chr "79P0060" "79P0061" "81P0725" "81P0729" ...
## $ pedon_id : chr "78IN177007" "78IN177008" "S1981IN011002" "S1981IN011006" ...
## $ pscs_top : int 15 30 20 20 25 15 NA NA NA NA ...
## $ pscs_bottom : int 65 80 70 70 75 65 NA NA NA NA ...
## $ pedon_completeness_index: int 1 1 7 7 6 7 6 6 6 6 ...
## $ state : chr "Indiana" "Indiana" "Indiana" "Indiana" ...
## $ mlra : chr "111A" "111A" "111A" "111A" ...
## $ geomposhill : chr "" "" "" "" ...
## $ geomposmntn : logi NA NA NA NA NA NA ...
## $ geompostrce : logi NA NA NA NA NA NA ...
## $ geomposflats : logi NA NA NA NA NA NA ...
## $ hillslopeprof : chr "" "" "" "" ...
## $ geomslopeseg : chr "" "" "" "" ...
## $ pmgroupname : chr "" "" "" "" ...
## $ drainagecl : chr "" "" "Well drained" "Well drained" ...
## $ taxonname : chr "Miami" "Miami" "Miami" "Miami" ...
## $ taxorder : chr "Alfisols" "Alfisols" "Alfisols" "Alfisols" ...
## $ taxsuborder : chr "Udalfs" "Udalfs" "Udalfs" "Udalfs" ...
## $ taxgrtgroup : chr "Hapludalfs" "Hapludalfs" "Hapludalfs" "Hapludalfs" ...
## $ taxsubgrp : chr "Typic Hapludalfs" "Typic Hapludalfs" "Oxyaquic Hapludalfs" "Oxyaquic Hapludalfs" ...
## $ taxpartsize : chr "fine-loamy" "fine-loamy" "fine-loamy" "fine-loamy" ...
## $ taxpartsizemod : chr "" "" "" "" ...
## $ taxceactcl : chr "" "" "active" "active" ...
## $ taxreaction : chr "" "not used" "" "" ...
## $ taxtempcl : chr "mesic" "" "mesic" "mesic" ...
## $ taxmoistscl : chr "" "" "Oxyaquic" "Oxyaquic" ...
## $ taxtempregime : chr "mesic" "mesic" "mesic" "mesic" ...
## $ soiltaxedition : chr "" "" "tenth edition" "tenth edition" ...
## $ osdtypelocflag : int 0 0 0 0 0 0 0 0 0 0 ...
## $ x : num -85 -85 -86.3 -86.3 -86.3 ...
## $ y : num 39.8 39.8 40 40 40 ...
## $ geometry :sfc_POINT of length 56; first list element: 'XY' num -85 39.8
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr "pedon_key" "pedlabsampnum" "pedon_id" "pscs_top" ...
# convert
# miami_sf <- st_as_sf(miami_sp)
# miami_sp <- as(miami_sf, "Spatial")
library(maptools)
st <- map("state", fill = TRUE, plot = FALSE)
# convert to sp object
st_sp <- map2SpatialPolygons(st, IDs = st$names)
proj4string(st_sp) <- CRS("+init=epsg:4326")
st_sp$state <- st$names
# loads data with ggplot2 package
st <- ggplot2::map_data("state")
library(broom)
st_tidy <- tidy(st_sp, region = "state")
IN_tidy <- tidy(as(IN, "Spatial"), region = "MUSYM")
ggplot2 plots data frames, but can also use sf objects (which are a special case of data frames). It can plot rasters, but only if they are converted to a data frame.
library(ggplot2)
# Lines
ggplot() +
geom_point(data = miami, aes(x = lon, y = lat)) +
geom_path(data = st, aes(x = long, y = lat, group = group)) +
xlim(range(miami$lon)) +
ylim(range(miami$lat)) +
ggtitle("Location of Miami Lab Pedons")
# Polygons
ggplot() +
geom_polygon(data = st_tidy, aes(x = long, y = lat, group = group, fill = id)) +
coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
# remove legend
guides(fill = FALSE)
# sf objects
# Polygons
ggplot() +
geom_sf(data = cnty, aes(fill = statefp, lty = NA)) +
geom_sf(data = miami_sf) +
coord_sf(crs = "+init=epsg:5070") +
guides(fill = FALSE)
# Facets
test <- subset(cnty_corn, year %in% 2012:2017)
ggplot() +
geom_sf(data = test, aes(fill = Value, lty = NA)) +
scale_fill_viridis_c(na.value = "transparent") +
facet_wrap(~ year) +
geom_path(data = st, aes(x = long, y = lat, group = group)) +
ggtitle(corn_yield$short_desc[1])
ggmap expands ggplot2 to download and plot base maps.
library(ggmap)
# build bound box and get base map via ggmap
bb <- sf::st_bbox(IN)
bb <- make_bbox(lon = bb[c(3, 1)], lat = bb[c(2, 4)])
gmap <- get_map(bb, maptype = "terrain", source = "osm")
# Lines
ggmap(gmap) +
geom_path(data = IN_tidy, aes(x = long, y = lat, group = group))
# geom_sf(data = IN, fill = NA, inherit.aes = FALSE) +
# guides(fill = FALSE)
# geom_sf() doesn't work with ggmp, their is a systematic shift https://github.com/r-spatial/sf/issues/336
The “t” in tmap stands for thematic, but tmap can also plot rasters natively. tmap’s syntax is very similar to ggplot2, but with a few twists.
library(tmap)
tm_shape(IN) + tm_polygons("MUSYM", border.col = NULL) +
tm_shape(cnty) + tm_borders() +
tm_shape(miami_sf) + tm_dots() +
tm_legend(legend.outside = TRUE)
# interactive web map
tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(IN) + tm_borders()
library(mapview)
cols <- RColorBrewer::brewer.pal(50, "Paired")
test <- mapview(IN, zcol = "MUSYM", lwd = 0, col.regions = cols) +
mapview(cnty, type = "l")
test
# export to html
mapshot(test, url = "C:/workspace2/test.html", selfcontained = FALSE)
library(leaflet)
test <- leaflet() %>%
addProviderTiles("Esri.WorldImagery", group = "Imagery") %>%
addPolygons(data = IN, fill = FALSE, color = "black", weight = 2)
test
# export to html
htmlwidgets::saveWidget(test, file = "C:/workspace2/test.html", selfcontained = FALSE)
Healy, K., 2018. Data Visualization: a practical introduction. Princeton University Press. http://socviz.co/
Gimond, M., 2019. Intro to GIS and Spatial Analysis. https://mgimond.github.io/Spatial/
Hijmans, R.J., 2019. Spatial Data Science with R. https://rspatial.org/
Lovelace, R., J. Nowosad, and J. Muenchow, 2019. Geocomputation with R. CRC Press. https://bookdown.org/robinlovelace/geocompr/
Pebesma, E., and R. Bivand, 2019. Spatial Data Science. https://keen-swartz-3146c4.netlify.com/